#include "SinglePhasePyrolysisADMaterial.h"
#include<iostream>
#include "libmesh/quadrature.h"
using std::cout;
using std::endl;
registerMooseObject("TrilobitaApp", SinglePhasePyrolysisADMaterial);

template<>
InputParameters validParams<SinglePhasePyrolysisADMaterial>()
{
	auto params = ADMaterial::validParams();
	params.addParam<std::string>("rho_name", "Material property density name");
	//add regular thermal properties
	params.addRequiredParam<Real>("rhov_i", "density of original material");
	params.addRequiredParam<Real>("rhoc_i", "density of fully charred material");
	//add special parameters to compute pyrolysis Arrhenius relation
	params.addRequiredParam<Real>("T_begin", "temperature at which the reaction begins");
	params.addRequiredParam<Real>("precofficient", "pre-exponential factor");
	params.addRequiredParam<Real>("ER", "E/R");
	params.addRequiredParam<Real>("reaction_level", "Reaction level of pyrolysis,normally 1,2 or 3");
	//add coupled variables
	params.addRequiredCoupledVar("temperature", "Coupled Temperature");
	params.addRequiredCoupledVar("rhoi", "Coupled density");
	return params;
}

ADReal SinglePhasePyrolysisADMaterial::Source(ADReal rho, ADReal T, Real T_begin, Real rhoc, Real rhov, Real precoff, Real ER, Real reaction_level)
{
	if (T > T_begin)
	{
		if ((rho - rhoc) / (rhov - rhoc) < 0.01)
		{
			return 0;
		}
		else
		{
			ADReal tmp = (rho - rhoc) / rhov;
			return precoff * exp(-ER / T) * rhov * pow(tmp, reaction_level);
		}
	}
	else
	{
		return 0;
	}
}
SinglePhasePyrolysisADMaterial::SinglePhasePyrolysisADMaterial(const InputParameters& parameters) :
	ADMaterial(parameters),
	_rho_name(isParamValid("rho_name") ? getParam<std::string>("rho_name") : ""),
	_pyrolysis_source(declareADProperty<Real>("pyrolysis_source" + _rho_name)),
	_T_value(adCoupledValue("temperature")),
	_rho_value(adCoupledValue("rhoi"))
{
	_rho_i_v_value = getParam<Real>("rhov_i");
	_rho_i_c_value = getParam<Real>("rhoc_i");
	_T_begin_value = getParam<Real>("T_begin");
	_precoff_value = getParam<Real>("precofficient");
	_ER_value = getParam<Real>("ER");
	_reaction_level_value = getParam<Real>("reaction_level");
}

void SinglePhasePyrolysisADMaterial::computeProperties()
{
	for (unsigned int qp(0); qp < _qrule->n_points(); ++qp)
	{
		_pyrolysis_source[qp] = Source(_rho_value[qp], _T_value[qp], _T_begin_value, _rho_i_c_value, _rho_i_v_value, _precoff_value, _ER_value, _reaction_level_value);
	}

}


